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I. INTRODUCTION 

One of the most important goals of the experimental program of high energy nuclear collisions is to determine the 
transport and thermodynamical properties of QCD as a function of temperature T and baryon chemical potential /i. 
During the past few decades, the experimental data measured at the Relativistic Heavy Ion Collider (RHIC) at the 
Brookhaven National Laboratory and the Large Hadron Collider (LHC) at CERN in Geneva, has shown collective 
behavior of the QCD matter created after the collision of heavy nuclei at high energies HHZI- The low momentum 
region of the transverse hadron spectra and the two particle correlation functions are well described by relativistic 
viscous fluid dynamics with a very small value of the shear viscosity over entropy ratio [/] These results have been 
taken as evidence for the production of an almost perfect liquid, a strongly coupled quark gluon plasma. 

The hydrodynamic modeling of heavy ion collisions solves on an event-by-event basis the relativistic fluid equations 
corresponding to energy-momentum conservation laws together with the so called constitutive relations for the shear 
viscous tensor and bulk pressure. Within this approach, little attention has been paid to the possible role of the 
baryon density n and/or baryon chemical potential fx. At high energies, this is justified because n and /r are very 
small, at least in the midrapidity region. However, interesting physics could be probed by investigating event-by-event 
fluctuations in the local baryon number density. 

Baryon number fluctuations have been mainly discussed in the context of heavy ion collisions at lower energy where 
larger values of fx can be realized. Interesting features of the QCD phase diagram can emerge there m- Different 
effective models have predicted the existence of a first-order phase boundary that separates hadronic matter from 
the quark gluon plasma at larger values of the baryon chemical potential. This boundary comes to an end at some 
critical values of the temperature T c and baryon chemical potential /r c . Right now there is no conclusive evidence for 
the location of a critical point in the T — /i plane from lattice QCD calculations at finite baryon density El- 

On the other hand, in heavy ion collisions it has been proposed to study second and higher order cumulants 
of particle multiplicity distributions as a function of the center of mass energy HMD- From thermodynamic 
considerations, it is expected that these moments scale with the correlation length which is expected to become large 
near the QCD critical point PEMH1H3 M EH- Possible signs of the critical point have been measured at RHIC 
but at present these do not provide a conclusive evidence [3314361 . If the expanding fireball of nuclear matter passes 
through a critical region (close to a critical point), one can extract information about the equation of state and the 
critical behavior of transport coefficients from the particle spectra formed at the freeze-out surface. It is important 
to determine whether the possible signatures of the critical point can survive the entire evolution of the expanding 
fireball. 

In the fluid dynamic framework, different aspects of the evolution of the fireball can change the pattern expected 
from purely thermodynamic considerations. Thermodynamic fluctuations are in principle part of a fluid dynamic 
description, at least in an extended sense where one accounts also for noise. Fluctuations evolve in time and space 
during the expansion of the fireball and thus, these are indeed effected by the equation of state and specially the 
transport coefficients such as the viscosities and conductivities. Close to equilibrium, there is also a deep theoretical 
connection between thermodynamic fluctuations in fluid dynamic fields and dissipative transport properties as stated 
by the fluctuation-dissipation theorem; see, e. g., Refs. [371439) . In the vicinity of the critical point, heat conductivity 
k as well as the shear and bulk viscosities rj and £ show critical behavior [4D14i5] . 

Besides genuine thermodynamic fluctuations (or noise), there is another possible source of fluctuations in the fluid 
dynamic approach to heavy ion collisions. These are the fluctuations already present in the initial state when the 
fluid dynamic treatment becomes valid. Their origin can be either the substructure of the colliding nuclei or the far- 
from-equilibrium dynamics preceding a fluid dynamic regime. This kind of initial state perturbation is particularly 
important for energy and/or entropy density. Fluctuations in the geometric distribution of nucleons within a nucleus 
lead to initial density perturbations which - after a fluid dynamical evolution - determine the spectrum of harmonic 
flow coefficients and, for example the form of the two-particle correlation function (“the ridge”) in heavy ion collisions. 

In a very similar way to fluctuations in the initial energy density, one can also expect, for example from a Glauber- 
type description of the initial state, initial fluctuations in the baryon number density. Indeed, baryon number density 
carried by protons and neutrons is presumably not distributed homogeneously within a nucleus and fluctuates locally 
and from event to event. In addition, the baryons and anti-baryons produced by pair production directly after the 
collision are subject to some local and event-by-event fluctuations [44] . 

In order to discriminate the effects associated to the thermodynamic fluctuations from the initial state fluctuations, 
it is necessary to understand their space-time evolution. In the present work we will concentrate mainly on the dynam¬ 
ics of initial state fluctuations although parts of our formalism are relevant also for the evolution of thermodynamic 


1 See Refs. mill for the recent developments in relativistic hydrodynamics. 
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fluctuations. Initial state fluctuations are interesting on their own. For instance, the evolution of the initial pertur¬ 
bations of energy density depends on the viscosities, in particular shear viscosity. In a similar way, the evolution of 
baryon number density depends on heat conductivity (in the Landau frame one may see heat conductivity equivalently 
as baryon number diffusion). If one has a theoretical understanding of initial state perturbations in baryon number 
density and their fluid dynamic evolution, it is possible to study their consequences for particle spectra at freeze-out. 
Provided possible signals are large enough to be seen within the constraints set by finite statistics, there could be a 
possibility to constrain the heat conductivity of the quark gluon plasma from experimental data. This would be very 
interesting for not only low energy collision experiments which aim at exploring the QCD phase diagram, but also at 
RHIC and LHC energies where baryon number diffusion could be another characteristic of the quark-gluon plasma. 

As a first step in this direction we study here the fluid dynamic propagation of local and event-by-event fluctuations 
of the baryon number density, flow velocity and energy density. These fluctuations propagate on top of a hydrodynam- 
ical background which for simplicity, we consider to be described by Bjorken’s model [45| (which includes finite baryon 
number density). In order to study the fluid dynamic propagation of perturbations we use a background-fluctuation 
splitting and a Bessel-Fourier decomposition for the fluctuating part of the fluid fields pRI]iTl| . We derive the evolution 
equations of the linear fluctuations and solve them for different initial conditions, values of the transport coefficients 
and equation of state. 

This work is organized as follows. In Sec. [IT] we review briefly the theory of relativistic fluid dynamics at finite chem¬ 
ical potential putting emphasis on the role of the equation of state and current estimates of the transport coefficients 
in the strong and weakly coupling regimes. The main features of the temporal evolution of the background fields are 
discussed in Sec. |III[ In Sec. ||V] we formulate the theory of linear perturbations on top of this evolving background 
and discuss numerical solutions. In Sec. [V] we draw some conclusions for a potential experimental observable, the 
correlation function of net baryon number as a function of azimuthal angles and rapidity. General conclusions are 
presented in Sec. VI Some technical details of our calculations are presented in Appendixes [A] and [B] respectively. 


II. RELATIVISTIC FLUID DYNAMICS WITH A GLOBALLY CONSERVED CHARGE 

We consider a relativistic fluid with one globally conserved quantum number current (baryonic number for our 
purposes). The energy-momentum tensor and number current are 


= e// + {p + 7T b ulk)A /iI/ + 7T^ , 


Here, e is the energy density, u^ is the fluid velocity, is the shear stress tensor, 7Tb u ik is the bulk viscous pressure, 
n is the particle density and zF is the particle diffusion current. We choose the signature of the metric to be 
(—, +, +, +) and the projector orthogonal to the fluid velocity is 

A'“' = f+tiV. (2) 


The fluid velocity is normalized to = — 1. We work in the Landau frame where the fluid velocity is chosen such 
that = —eu u . The shear stress tensor is transverse to the fluid velocity, 

u^ v = 0. (3) 

The shear stress tensor is also symmetric and traceless. The particle number density is defined by n = —u^N 11 such 
that the diffusion current is orthogonal to the fluid velocity, = 0. 

It is clear that an arbitrary (symmetric) energy-momentum tensor (with a time-like eigenvector) and current 
N M can be written in the above form. The decomposition becomes unique by requiring that the pressure p is related 
to the energy density e and the baryon density n by the same relation as in thermodynamic equilibrium, i.e. by an 
equation of state p = p(e, n) . 

The evolution of the energy momentum-tensor and the particle current are constrained by the conservation equations 


= 0 , 

V m A^ = 0, 


( 4 ) 


where denotes the covariant derivative. In this general form the conservation equations hold also in curved space- 
time but we are interested here in curvilinear systems defined in Minkowski space without taking into consideration 
the gravitational field. From Eqs. ([I]) and ([4]) one obtains the evolution equations for the energy density, fluid velocity 
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and particle density 


De + (e + p + 7r bulk )V M it M + = 0, 

(e + p + 7Tbulk) DvT + <9^ (p + 7Tbulk) + A 1 " a V/i7T^ a = 0, 

-D77 + nV^u^ + V m i/ m = 0 . 


( 5 ) 


Here we have introduced the comoving derivative defined as D = it M V M . 

To close the evolution equations ([5| one needs expressions for 7Tb u ik, 7 r M1/ and i/ M . Within the formalism of fluid 
dynamics one writes these objects as a derivative expansion in terms of the fluid velocity u^ and thermodynamic 
variables e, n. In the present work we concentrate for simplicity on the first order of this expansion. One should keep 
in mind that terms of second order are expected to improve the results quantitatively and are in general needed for 
an acceptable causal structure and linear stability O ‘55j. 

The constitutive relation for the shear stress is 


— 


= —2r] cr M " = — 2?? 


-A^A"' 8 + -A^A"“ — Ia^A^ 

2 2 3 


^ a^/3 5 


( 6 ) 


where 77 is the shear viscosity transport coefficient. The bulk viscous pressure is obtained from the following expression 


T^bulk — £# — CV fj,U^ , 

where £ is the bulk viscosity and 6 is the expansion scalar. Finally, the particle diffusion current is 

2 r rn I 2 


v = —K 


nT 

e+p 


l = —K 


nT 
e + p 


A' 


% (£) • 


(7) 


( 8 ) 


where k is the heat conductivity. In the last equation we have introduced the chemical potential p, which is conjugate 
to the baryon density n, and the temperature T. 

In summary, the hydrodynamic equations at this stage involve the fluid velocity u M (with three independent com¬ 
ponents), the energy density e, pressure p, baryon density n, baryon chemical potential p, temperature T as well 
as the shear viscosity 77 , bulk viscosity £ and the thermal conductivity k. Only two thermodynamic variables are 
independent and they also determine the transport properties 77 , £ and k. In a non-equilibrium situation only energy 
density e = and baryon number density n = are directly related to the physical energy-momentum 

tensor T and number current . All other thermodynamic variables are defined indirectly via their relation to e 
and n in thermal equilibrium. 

For the practical calculations one is in principle free to use any set of independent thermodynamic variables. The 
form of Eqs. ([5| suggests the use of the energy density e and baryon density n. However, because most microscopic 
calculations are done in the grand canonical ensemble, the thermodynamic equation of state and the transport 
coefficients are usually obtained as a function of the temperature T and chemical potential p, for example p = p(T, p). 
Thus, it can be advantageous to use T and p as independent variables in fluid dynamics, as well. This avoids the 
inversion of functions which can be numerically difficult. One should keep in mind that T and p in a non-equilibrium 
situation are defined via their relation to e and n. Eq. ([ 5 ]) can be transformed using thermodynamic relations compiled 
in Appendix |Aj The evolution equation for energy density becomes 


d 2 p d 2 p 
3T 2 ^dTdp 


DT 


T 


d 2 p d 2 p 
dTdp +M <V 


Dp + (e + p) 6 - 277 a a p(J a P - £ 0 2 = 0. 


(9) 


where we have now used the constitutive relations © and Q . The evolution equation for the fluid velocity is now of 
the form 

(e + p) Du v + A va {s d a T + n d a p ) - \\Vp (2 77 u Q/3 + £ A a ? V 7 u 7 ) = 0, (10) 

and finally, the particle number conservation law becomes 

J^ OT+ g^ + „ fl + v ^ = 0 . ( 11 ) 

Note that Eqs. © and form a linear system of equations that can be solved for DT = u a d a T and Dp = u a d a p 
as long as 


d 2 p d 2 p f d 2 p 
dT 2 d^~ 


7^ 0 . 


dTdp 


( 12 ) 
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To solve the fluid dynamic equations we will use a background-fluctuation splitting. To this end we write the fluid 
dynamic fields as 


+ 5u^, e = e + 5e, (13) 

and similar for the other fields. We are interested in perturbations Se etc. that are small enough so that only 
linear terms in the evolution equations need to be kept. The background fields u^, e etc. satisfy the fluid dynamic 
equations <[5j) while the perturbations satisfy linear equations that depend on the background solution. We derive these 
linearized equations for arbitrary background fields in Appendix [Bj The structure of the linearized equations permits 
us to simply use de, Sn and three independent components of the fluid velocity as variables (the fourth component 
of the fluid velocity follows from the constraint u^Su^ = 0). However, all the background-dependent thermodynamic 
quantities can be expressed in terms of T and p. Useful thermodynamic relations for this purpose are compiled in 
Appendix [A] 

In the rest of this section we briefly discuss some simple parametrization of the thermodynamic equations of state 
p(T,p) and transport properties rj(T,p), £(T, p), n{T,p). We emphasize that our formalism can be used for an 
arbitrary form of these functions once these have been determined from a particular microscopic description. 


A. Equation of state 


The fluid hydrodynamical equations require an equation of state (EOS) p(T, p) as an input. In principle, the 
equation of state can be calculated from the inherent quantum field theory associated to a particular system but this 
is a formidable task. In recent years there have been important advances to determine analytically and numerically 
the thermodynamical properties of QCD at high temperatures and chemical potential by considering effective thermal 
field theories |56W62| while in the low temperature and chemical potential regimes one expects that a non-interacting 
hadron resonance gas provides a reasonably good approximation [_63]. 

At intermediate temperatures, non-perturbative methods are needed to describe the transition which separates the 
hadronic, confined phase and the quark-gluon plasma (QGP) phase. While several studies of lattice QCD simulations 
are available at the moment at vanishing chemical potential p = a,t p > 0 lattice simulations are not possible 
due to the sign problem. However, different alternatives have been studied in order to circumvent this problem such 
as reweighting |64j . Taylor-expansion in p [651172) . analytic continuation from imaginary p [73], the density of states 
method, or using the canonical ensemble. Of course, each of these methods have their advantages and disadvantages. 

One of the main goals in the analysis of fluid dynamic fluctuations and their propagation is to provide a phenomeno¬ 
logical determination of the equation of state (or at least of some of its properties). In the derivation of the evolution 
equations for the background and fluctuating fields we shall keep the equation of state p(T , p) unspecified as far as 
possible in analytic expressions. For some numerical calculations and illustrations we use the simplest possible case, 
a non-interacting gas of Np massless quarks that come in Nq colors and Nq — 1 gluons, 


where we use the abbreviations 


P( T , P) = ^ a i tA + \ a i T 2 p 2 + ^j a 3 p 4 , 


a i 

d2 

a-3 



1 + -NcNp 


2 N c N f 
27 

2 N c N f 

8l7T 2 


(14) 


(15) 


The baryon chemical potential p measures the net baryon density of the system. In our convention, quarks carry 
baryon number charge 1/3 and anti-quarks —1/3. 

Corrections to the ideal EOS arise as a consequence of interactions and the breaking of conformal invariance by 
dimensional transmutation and non-zero quark masses. They are most important at low temperatures. We follow 
here the Wuppertal collaboration which has parametrized the QCD equation of state for finite chemical potential in 


2 For a recent review of the lattice QCD studies we refer to the reader to Ref. m- 
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Transport coefficient 

Weakly-coupled QCD 

Strongly-coupled theories 

V 

T 3 

S (7» 

C 9 4 log(l/ g) 

47r 

C 

15 V (T)(l-c 2 a (T)) 2 

2r 1 (T,g)(\-c 2 s {T,g)) 

K, 

~ /A /g* for g » T 
~ T 4 /( 5 y )for g « T 

8n 2 ^ V (T,g) 


TABLE I. Estimated values of the shear viscosity and different parametrizations for the bulk viscosities and heat conductivity 
for weakly-coupled QCD 17514791 and strongly coupled theories with holographic duals [8044821 . See text for discussion. 


terms of a Taylor expansion ED- The leading order expression for the trace anomaly or QCD interaction measure 
I(T,n) = e(T,g) - 3 p(T,g) is 


/(?» I(T, 0) g 2 d X2 (T) 

T 4 T 4 + 2 T dT 

where J(T, 0) is the interaction measured at g = 0 and X 2 {T) is the leading-order Taylor coefficient. 
J(T, 0) and X2 (T), can be parametrized analytically as [7TJ 75] 


(16) 

Both terms, 


I{T,0) „_ h , /t _ h7/t i \ u , / 0 (tanh(/i t + f 2 ) + 1) 

T 4 " 6 ' r° + 1+ 91 t + g 2 t 2 

X2 (T) = e-^A-^/t 2 h (tanh(/ 4 t + f 5 ) + 1) 


(17a) 

(17b) 


where t = T/( 0.2 GeV). For Nf = 2 + 1 flavors of quarks with physical masses and finite baryon chemical potential 
g the parameters in Eq. flT7| are h 0 = 0.1396, hi = 0.1800, h 2 = 0.0350, /o = 2.76, /i = 6.79, f 2 = 5.29, g\ = 0.47, 
g 2 = 1.04, h 3 = —0.5022, h 4 = 0.5950, f 3 = 0.0940, / 4 = 6.3290 and = —4.8303 [71 ( [ 71 ]. The pressure at finite g 
is given by 


P(T, g) = p{T , 0) 1 g^_ 

r j'4 y ’ 4 ' 2 T 2 ^ ~ ' 

At g = 0 the relation between the pressure and the trace anomaly (|17a) is 


P(T, 0) 

jn 4 


dT 


I(T', 0) 

Y 1 / 5 


(18) 


(19) 


All other thermodynamic quantities can be derived from p(T , g) using the standard relations (compiled in appendix 
0. The equation of state ( |18[ > with the above parametrization is valid for small chemical potentials g/T < 3 in the 
temperature window 0 < T < 400 MeV. We will use Eq. (18) to study the influence of the EOS for the dynamics of 
the background fluid dynamic fields. 


B. Transport coefficients 

In addition to the thermodynamic equation of state, the fluid dynamical description needs as an input transport 
coefficients. These can either be determined experimentally, or, if a microscopic underlying theory is known, they can 
at least in principle be calculated as a function of the thermodynamic variables via Kubo relations. In this section 
we briefly summarize the current theoretical knowledge for the shear and bulk viscosities and thermal conductivity of 
QCD and related theories, both in weakly and strongly coupled regime^] 


C. Weak coupling regime 

When the interaction strength is small, effective thermal held theory methods allows us to calculate the transport 
coefficients. For weakly coupled QCD in the high temperature and vanishing chemical potential regime, the leading 


3 A more detailed discussion of the properties of the transport coefficients discussed in this work can be found in Ref. m- 
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logarithmic result for the shear viscosity is (75H77] 


p(T) = k 


g 4 iog(i/ff) ’ 


( 20 ) 


where g is the strong coupling constant. In the previous expression k is a constant that depends on the number of 
fermions species [75HZ3- Arnold et al. showed that at leading log accuracy and for high temperatures with vanishing 
chemical potential there is an approximate scaling between the shear (p) and bulk (£) viscosities for weakly coupled 

QCD [77] 

C(T)«15»7(T)Q-c2(r)) , (21) 

where c 2 = dp/de is the speed of sound. A similar expression was first derived by Weinberg for a gas of photons [84] . 
To date there is no complete leading logarithmic calculation of the heat conductivity k(T, p) and so far only two 
estimates of k(T, p) have been provided in the literature for different kinematic regions of the T — p plane [75] [75] 


k(T, p) 


F[T , m D ) p 2 /g 4 , for p > T, 
C T 4 /(g 4 p 2 ), for /r « T. 


( 22 ) 


where F(T, mo) is a function that depends on the temperature and the Debye screening mass mo (see Ref. [751 for 
details). In the case of small chemical potential, the proportionality constant C depends on the number of flavors and 
the gauge group m In the limit where p —> 0 the heat conductivity n ~ p 2 is divergent. However, the particle 
diffusion current ([8]) remains finite CZ2I- In the context of relativistic kinetic theory, some general expressions for 
the transport coefficients with constant cross section or within the relaxation time approximation have been derived 
recently [55M50] . However, these calculations do not take into account the quantum screening effects of the QCD 
plasma. 

Despite relatively large uncertainties, experimental results indicate that the value of the shear viscosity over the 
entropy ratio p/s is smaller than the one calculated from weakly coupled QCD (20) [8)[9]. For the case of the bulk 
viscosity the situation is less clear: the uncertainties in its experimental determination are even larger (see Ref. m 
and references therein). In addition, there are no experimental constraints for the value of heat conductivity in high 
energy-nuclear collisions so far. 


D. Strong coupling regime 


From the previous discussion it is clear that at this moment perturbative QCD calculations of the transport 
coefficients are not completely under control for all the possible physical values of the temperature and chemical 
potential. On the other side, there are certain classes of strongly interacting theories where transport coefficients can 
be determined for almost all values of T and p. These are field theories with known gravitational duals where the 
computations can be done via the anti-de Sitter/conformal field theory (AdS/CFT) correspondence. Despite the fact 
that those theories are not equivalent to QCD, they share some qualitative aspects with it and thus, these theories 
might provide some guidance in the regimes where pQCD calculations are not reliable. We take here a pragmatical 
approach and consider the estimates of the transport coefficients based on holographic calculations as toy models 
which allow us to study the propagation of perturbations in fluid dynamic fields. For large t’Hooft coupling and for 
Af= 4 SYM theory, holographic methods give the well known result [50] , 


v(T,p) = J_ 

s(T, p) 4tt ' 


(23) 


This result holds also for any holographic theory at sufficiently large coupling and number of colors as long as the 
theory is spatially isotropic. This relation for p/s holds even in the presence of non-zero chemical potential [81] . 
Initially this result was conjectured to be an universal lower bound but today there is evidence showing that this 
relation does not hold in general [52H55] . Incidentally, the value of the shear viscosity extracted from experiments in 
high energy nuclear collisions is closer to the one predicted for strongly coupled theories (23) than the one calculated 

(see Ref. [8] for a recent review). 


in weakly coupled QCD ([20]) 

The shear viscosity has also been calculated for pure Yang-Mills theory using lattice gauge theory for specific values 
of temperature [ 1001(101] . The estimated values for p/s are somewhat above the AdS/CFT values. Similarly, p/s as 
a function of temperature for vanishing baryon chemical potential has also been estimated for Yang-Mills theory as 








well as QCD by using diagrammatic functional relations and gluon spectral functions obtained by numerical analytic 
continuation from Euclidean quantum field theory [102111031 . The minimal value for QCD was found to be rj/s ss 0.17 
at temperature T ~ 1.3 T c . 

For holographic theories that deviate from conformal behavior the bulk viscosity has also been calculated 


C(T, f j) = 2r i (T, li ) ( ^-c 2 (T,/x) 


(24) 


As in the case of the shear viscosity value (23) this relation holds for certain theories with finite chemical potential [105] 
but it is not an universal bound m ■ By comparing the scalings between ( and 77 , Eqs. ( ]24[ ) and ( |2l| ), one observes 
that they differ in the strong and weak coupling regime. This mismatch between both parametrizations is currently 
not understood. In the case of the thermal conductivity k, the calculations for strongly coupled plasmas with finite 
chemical potential give the following result [5T| 


k(T,/x) = 8tt 2 -^7?(T,/x) , 


(25) 


which is an analog of the Wiedemann-Franz law [107| 4 5 [ As in the weakly coupled case (22), the heat conductivity is 


divergent ~ g 2 while the particle diffusion current (81) is finite. Recently the temperature-dependence of the first 
and second order transport coefficients have been studied in a particular holographic model m- 

We summarize the discussion presented in this section in Table [IJ where we show the estimates of the transport 
coefficients in both strong and weak coupling regimes. Mainly for reasons of simplicity, we shall concentrate here 
on the parametrizations of the transport coefficients in the strong coupling regime Eqs. (23), (241 and (p5|) for our 


numerical calculatons. Another advantage of using the parametrization of strongly coupled theories is that both 
transport coefficients, the bulk viscosity £ and the heat conductivity k, are proportional to the shear viscosity 77 and 
thus, one can not only study the effect of the dissipative corrections but also one can investigate the ‘weak’ and 
‘strong’ regimes by varying the values of 77 /s. We keep the functions g(T,n), £(T, fi), k(T, fi) unspecified as far as 
possible in our analytic calculations. 


III. BJORKEN BOOST INVARIANT SOLUTION 


In this section we study the solutions of the fluid dynamical equations for a quark-gluon plasma undergoing boost 
invariant longitudinal expansion. We assume translational and rotational symmetry in the transverse plane and arrive 
at a simple model for the early stages of a heavy ion collision first studied by Bjorken [45]. Our analysis is extended 
to the case where there is a non-vanishing baryon number density. The relatively simple homogeneous solutions will 
also serve as a background for a more elaborate discussion of perturbations around it in Sec. El 

It is convenient to change from Cartesian coordinates x M = (t,x 1 ,^ 2 , ^ 3 ) to the Milne coordinates (r, r, </>, 77 ) where 
T = — x\ is the longitudinal proper time, 77 = arctanh(x 3 /t) is the longitudinal (space) rapidity and r and </> are 

the usual polar coordinates in the transverse plane. The metric in the Milne coordinates is g^ = diag(—1, l,r 2 ,r 2 ). 
The main advantage of using these coordinate systems is that the symmetries of the Bjorken solution are explicitly 
manifest. Specifically, the symmetry group ISO(2)®SO(l, 1 )®Z 2 consists of translations and rotations in the transverse 
plane, longitudinal boosts 77 —> 77 + A 77 and reflections 77 —> -77 m- The Bjorken flow velocity = (1, 0,0,0) is the 
only invariant unit vector and the symmetry also implies that all fluid dynamic fields depend only on the longitudinal 
proper time r |45|. 

From Eqs. (pH) one finds that the evolution equations for energy density and particle number density are 


d T e + (e + p) 


1 

r 


(1*7 + 0 -2 =0, 


d T n + n— = 0 . 
r 


(26) 


We have used here the Christoffel symbols of the Milne coordinate system. The non-vanishing ones are T(L = 
T= 1 /t, T t w) = r, T+ = T* r = 1/r, T^ = —r. The shear tensor defined in Eq. ([ 6 ]) becomes a^ = 


4 We pointed out to the reader that Eq. \2A\ was derived by means of the gauge/gravity duality in Ref. m for a specific model. Other 
non conformal field theories El where the duality holds provide some modifications to the parametrization given by Eq. ( |24| >. 

5 The relation ( |25| > was derived originally for a conformal holographic theory. However, this expression does not hold for non-conformal 
systems within the AdS/CFT correspondence Eon. 
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diag (0, — — 3 ^ 2 , 373 ) with a^cr^ 1 ' = The expansion scalar is 9 = \ and the projector orthogonal to the 

fluid velocity is = diag(0,1,1,1). The particle diffusion current ([ 8 ]) is a vector orthogonal to u^ and therefore 
vanishes exactly for the Bjorken flow. 

While the particle number density is simply diluted by the one-dimensional expansion, the evolution of energy 
density in (261 contains an additional loss term from the thermodynamic work done by the expansion and a gain term 
from shear and bulk viscous effects. After the variable change to T and /i eq. (26) becomes 

(1 _ 477/3-kA d 2 p 
sTt ) dp 2 


d T T + 


n d 2 p 
r dTdp 


d 2 p d 2 p f d 2 p \ ‘ 

WT 2 Tip 2 \dTd]lJ 


= 0, 


d T fi • 


n d 2 p 
r dT 2 




1 _ 47//3+C \ d 2 y 
1 sTt I dTdp 


(27) 


d 2 p d 2 p _ / d 2 p \ * 


= 0 . 


dT 2 dp, 2 ^ dTdp J 

We observe that the size of viscous corrections to an isentropic expansion is determined by the parameter 

477/3 + c 


7 = 


sTt 


(28) 


Formally, the gradient expansion underlying viscous fluid dynamics can be used for 7 < 1 . Note that for a given 
thermodynamic equation of state p(T,n) and viscosities ?y(T, p.), £(T, /i) one can solve the two coupled ordinary 
differential equations (27). 

In the remainder of this section we discuss as a simple illustrative example the equation of state of an ideal gas of 
massless quarks and gluons in Eq. (14). The evolution equations (271 for the temperature T and chemical potential 
p, become 


1 


d T T + —T - — 


3t 


yT\ |aia 2 T 4 + (|aia 3 + a|)T 2 /7 2 + a 2 a 3 p 4 


1 


~ 


aia 2 T 4 + (aia 3 — 3 a\)T 2 p? + a 2 a 3 7 i 4 
|aia 2 T 3 7r + 2a|T/i 3 
r J aia 2 T 4 + (aia 3 — 3a|)T 2 /i 2 + a 2 a 3 7i 4 


7 T 


= 0 , 


= 0 , 


(29) 


where the coefficients a 1 , a 2 and a 3 are given in Eq. (15). Note that we use conventions where fj, is the chemical 
potential for baryons, and the chemical potential for quarks is p q = /r/3. 

Let us first discuss some interesting limiting cases of Eqs. (29): 

1. Ideal fluid dynamic expansion. When shear and bulk viscosities vanish, 77 = £ = 0, the temperature and 
the chemical potential decouple from each other. This allows us to solve Eqs. (29) exactly, which gives 


T(t) = T(t 0 )(^) 1/3 
M t ) = m(to) (y) 7 


(30a) 

(30b) 


The scaling solution of the temperature is not modified by the presence of the chemical potential and it coincides 
with the well known result found by Bjorken [35] , 

2. Vanishing chemical potential. The point with fi = 0 corresponds to a (partial) fixed point of the evolution 
equations (29) with extended symmetry (baryon number parity). The evolution equation for temperature 
becomes 


d T T + 


T 
3 r 


(l-7) = 0. 


(31) 


where 7 is given ^ 
the previous equation is 


(28). For vanishing bulk viscosity, £ = 0, and constant ratio 77 /s, the exact solution to 


TTTt-fTT4j 


T(t)=T(t 0 ) (^) 


1/3 


1 + 


3r 0 T(r 0 ) s 


-(!-(— 


(?) 


2/3 


(32) 


Viscous corrections are relevant only at early times where velocity gradients are large while at late times these 
are suppressed and thus, T(r) ~ r -1 / 3 . 
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3. Small chemical potential. For [x/T -c 1 the dynamics of T is approximmately determined by Eq. (31) while 
the evolution equation for /i is 


d T n + ^ (1 + 2y) — 0. (33) 

The viscous effects (encoded in the parameter 7 ) have the tendency to accelerate the decrease of jx due to the 
expansion. This is in contrast to the temperature where viscosity has the opposite effect. To lowest order in 
77 /s, the solution of |33|) is 


K t ) = Mn>) (^) 


1/3 


1 - 


4 77 / _ n _ o \ 2 / 3 

3roT(r 0 ) s V ' t ) 


(34) 


4. Small temperature. For T/fj, <C 1 the evolution equation for the chemical potential is the one of eq. (30b) 
with a simple scaling solution. For the temperature we obtain to lowest order in T / fx 


drT+^( l-3 7 )=0, 


(35) 


which has a solution similar to Eq. ([32] ) when 77 /s and <//s have constant values. If one chooses T(r 0 ) = 0 as 
initial condition the solution to Eq. 1351) becomes 


T(t) = 


477+ 3 C 

2 s 


1 


1 


3 T t 


(36) 


Even if the temperature vanishes initially, the system is heated up due to shear and bulk dissipative effects. In 
contrast to [x = 0, vanishing temperature T = 0 does not correspond to a (partial) fixed point of the evolution. 


Let us now consider the evolution equations (|29j) in the general case where we find their solution numerically. In 
Fig. 0 we show the time evolution of the temperature (left panel) and chemical potential (right panel) for different 
constant values of 77 /s = 0 and 77 /s = 2 /( 47 r) (black and red lines respectively) and two different parametrizations of 
the equation of state: the ideal EOS (141 (solid lines) and the lattice-based EOS (fl 8 [ (dashed lines). The initial values 


at time r 0 = 0.5 fm/c are taken to be T{tq) = 0.4 GeV and fx(r 0 ) — 0.4 GeV. For the ideal EOS the coefficients an, 
02 and 03 are taken according to Eq. (15) with Nc = Np = 3. The bulk viscosity (24) vanishes exactly for the ideal 


EOS (14) but it becomes a function of the temperature and chemical potential for the lattice-based EOS (18). 


First we discuss the properties of the numerical solutions of Eqs. (27) for the ideal (and massless) EOS (14). For 
both variables, T and /x, the effect of viscous corrections are more relevant during the early stages of the expansion 
while at late times their effects are negligible as expected. In the left panel of Fig. [Tjwe see that the viscosity reduces 
the effect of the longitudinal expansion on the temperature. This is simply the expected heating by dissipative effects. 
At the final time 17 = 10 fm/c the temperature is larger by values of the order of 10% for 77 /s = 2 /( 47 r) compared 
to the ideal fluid expansion. For the chemical potential we find that the inclusion of dissipative corrections has the 
opposite effect, i.e. the chemical potentials decrease faster in the viscous case. This is clearly seen in the right panel 
of Fig. [l] when comparing the final values of the chemical potential 77 ( 77 ). The changes with respect to the ideal fluid 
expansion are also somewhat larger, of the order of 15% for 77 /s = 2/(47 t). 


When using the lattice-based EOS (18 1 we find that the numerical solutions of Eqs. (27) for T and /i are qualitatively 


similar to the ones obtained from the ideal EOS during the early stages of the evolution. As a function of time, the 
temperature is always decreasing and the dissipative corrections are larger at early times than at late times. The 
chemical potential decreases faster for larger values of the shear viscosity. For the lattice EOS, the changes induced 
by the dissipative corrections are on the order of 8 — 15%. 

Interestingly, the evolution of /t with time differs substantially between the two choices for the equation of state. 
In the right panel of Fig. [j] one observes that the decrease with time is much weaker for the lattice EOS than for the 
ideal EOS. At the freeze-out time 17 = 10 fm/c and for vanishing 77 /s, one has 77 ( 77 ) ss 0.29 GeV for the lattice EOS 
while ix(rf) ~ 0.12 GeV for the ideal EOS. The difference between those values increases slightly for finite values of 
77 /s. Moreover, at late times /r increases slowly (and somewhat more for larger values of 77 /s). Our numerical results 
show also that when using the lattice EOS the values of the temperature are somewhat larger than for the ideal EOS 
specially at late times. 

In Fig. [2] we show the Bjorken flow trajectories in the plane of chemical potential 77 and temperature T for the 
ideal EOS (|14| (left panel) and the lattice EOS (right panel). For the initial conditions we choose r 0 = 0.5 fm/c, 
T(r 0 ) = 0.4 GeV and different values of /x(to) = {0.05, 0.15, 0.25,0.35,0.45, 0.55} GeV. For both equations of state we 
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(a) (b) 



FIG. 1. (Color online) Log-log plot of the time evolution of (a) temperature (left panel) and (b) chemical potential (right panel) 
for the ideal EOS (141 (solid lines) and the lattice-based EOS (181 (dashed lines). We choose here rj/s = 0 (black line) and 
r//s = 2/(Atv) (red line). For the initial conditions we select T(tq) = fi(j o) = 0.4 GeV and tq = 0.5 fm/c. 



FIG. 2. (Color online) Flow trajectory of Bjorken expansions in the /r-T-plane, initialized for (a) ideal EOS ( |14[ ) (left panel) 
and (b) lattice-based EOS ( |18| ) (right panel). For the initial conditions we choose to = 0.5 fm/c, T{tq) = 0.4 GeV and different 
values of /t(to) = {0.05, 0.15, 0.25, 0.35, 0.45, 0.55} GeV. In both panels we compare the viscous effects by choosing tj/s = 2/(47r) 
(dashed lines) to the the case of vanishing viscosity, rf/s = 0 (solid lines). All lines end at fixed final time r/ = 10 fm/c. Note 
that we use conventions where fi is the chemical potential for baryons, the chemical potential for quarks is fi q = fi/3. 


vary the shear viscosity to entropy if/s — 2/(47 t) (dashed lines) and rf/s = 0 (solid lines). All trajectories end at fixed 
final time r = 10 fm/c. 

For the ideal EOS (left panel of Fig. [2]) we observe that the viscosity weakens the effect of the expansion on the 
temperature T while it does the opposite for the chemical potential fi and thus the trajectories end at larger values of 
T and smaller vales of fi for non-zero if/s. This is in agreement with the previous discussion of the temporal evolution 
of T and fi. For the lattice EOS (right panel of Fig. [2]) we observe similar trajectories for small initial values of fi(r 0 ). 
For larger values of ^(r 0 ), the trajectories start to bend towards larger values of fi while they continue to decrease 
towards lower values of T. This behavior is understood from the previous discussion, as well. 

In summary, the time-evolution of temperature and chemical potential for a Bjorken expansion is given by Eqs. (27) 
for an arbitrary EOS. The evolution of p, as a function of time is quite sensitive to the choice of the EOS. The effect 
of the viscosity is relatively small. This is actually expected for the homogeneous background while we expect more 
prominent dissipative effects for non-homogeneous perturbations around it. [^] We turn to those in the next section. 


6 The effect of shear viscosity is also sizable for the transverse expansion (radial flow) and for elliptic flow m- 
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IV. FLUCTUATIONS AROUND BJORKEN FLOW 


After having studied the solution of the hydrodynamic evolution equations with Bjorken boost invariance and 
transverse translational symmetries we study now the evolution of fluctuations or deviations from that solution. We 
will concentrate here on deviations that are small enough in magnitude to describe their evolution by linearized 
evolution equations. In other words, we write the fluid dynamic Helds as 


+ Su^, e = e + 6e, n = n + Sn, 


(37) 


where u 11 , e, n is the Bjorken-type solution discussed in the previous section. The linearized evolution equations for 
the perturbations SuSe, Sn are discussed for a generic background solution and arbitrary coordinate system in 
Appendix[B] If one specializes to the Bjorken background and the coordinate system (r, r, </>, rf), the independent fluid 
dynamic fields are in the first order formalism Se, Sn, Su r , Su 4 and Su v . (We take the background fluid velocity u M 
and the full fluid velocity + Su^ to be normalized, u^u^ = = —1, such that one has Su T = 0 at linear 

order in perturbations). Equation (B3) yields the following equation for the perturbation in energy density (each 
hydrodynamical fluctuating field depends on (r, r, <j>, rf) which we suppress for better readability) 


d T Se + 

+ 

+ 


1 

T 


-(if) -4(f) 

T \ de Jn T \ de J 


4 

3f2 



r \dn) e t 2 \dn) e 


4 

3 T 2 



Sn 




d r Su r + -5u r + dju 4 + dju 71 
r 


-rj drjSu 11 = 0 . 


(38) 


The thermodynamic derivatives like (dp/de) n , etc., are to be evaluated here on the background solution and similarly 
the transport coefficients and their derivatives. The evolution equation for the perturbation in baryon number density 
is 


d T Sn H —Sn 

T 


nT " 2 
e+p 


*(£) 


nT n2 


e + p 


nT " 2 


e + p 


d r Su r -\ — u r + d^Su 4 + d^Su 71 
r 

(+aP), (5 + l a - + yA + %) *> = o- 


(39) 


The derivative operator of second order that appears in the last two lines in front of Se and Sn, respectively, is the 
Laplace operator in the spatial coordinates r, </> and 77 . 

The fluid velocity in the radial direction is determined by the following evolution equation 


e+p - -C + p-p ) d T Su r + 
t 3r 


d T p - ^d T ( +^C+ ^ d T p + ^77 


Su r 


+ 

-c 

- v 


dp 

~de 


1 (d( 

de 


dr/ 

3 t V de 


d r Se - 


dp 

dn 


1 fd( 

dn 


dp 
3 t V dn 


d r Sn 


df + -d r —^ ) Su r + drd^Su 4 + drd^Su 71 


+ i-a r - A + _L 9 2 + §u r + (}_ drdtj) _ 5u 0 + 1 drdri§u . 


3 ? ,2 1 r 2 <t> ' r 2 “V 


= 0 , 


(40) 
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the one in the azimuthal direction by 


e +p - -C + p ) d T Su^ + 


3 T 

1 (d C 

~d~e 


d T p - -d T C, + drV + 


-c 


1 


dp 
3 r V de 
1 


12 dc/,Se + 


3 r 
dp 
dn 


3 r 2 

1 /<9C 




t \dn 


dp 
3 r V dn 


^2 d <f>Sn 


-^dj + —5^ 5u r + —d^Su 4, + —dJSu 1 


3 r 2 


3 r 3 


1 


5,9. + „ J d> Uu r + dl + -5, + —:dl + Ji Su 4 + 


3r 2 0 


r 2 V 


^d+dnSu" 


= 0 , 


and finally the fluid velocity component in the rapidity direction is governed by 


e +P - ) drSu v + 


2 1-1-4 4 

9,P+ -(e + p)- drC + - —9 r ?y — 


9p 

9^ 


-C 

- fj 


1 


I ^ 

r \9e 
1 


dp 
3r V 5e 
1 


~^d v 8e + 


9p 

9n 


~^d r d v + d v 8u r + —dj v 5u 4 + —dlSvP 


1 


3r 2 


3r 2 r 


1 


9 r 9„ + „ d r ,) Su r + —.djSu 4 + [di + -d r + Ji + „ Ji) Su ri 


3 T 2 


1. 


3 r 

1 3 

r V dn 


1 


3 T 2 
4 

3t V 9n 


<5^ 


9ry 


r 2 <t> — 3 t 2 ^ 


^ 2 dr,8n 


= 0 . 


(41) 


(42) 


Equations (381 - (42) are hyperbolic coupled linear differential equations for the variables 8e , Sn , 5u r , Su 4 and Su v : 
They contain only first order derivatives with respect to the time coordinate r but up to second order derivatives 
with respect to the spatial coordinates r, (j) and p. In the second order gradient expansion the equations would be 
elliptical but also contain more degrees of freedom and transport coefficients. 

In order to analyze the differential equations (38) - (42) it is convenient to use a Bessel-Fourier transformation. For 
the perturbation in energy density this reads 


Se{T,r,(/),p) = J dkk ^ J ^ Se(r, k, m, q) J m (kr), 


(43) 


with inverse relation 


POO 1 p2,7T P 

Se(T,k,m,q) = J drr— j d<j> J dp Se(r, r, (f>, p) e~ l( ' m4>+qv ' > J m (kr). 


Since Se(r,r, cj),p) £ M and J- m (kr) = (—1 ) m J m (/cr) one has 

8e*(r, k, to, q) = (—1 ) m <5e(r, k, —to, —q). 


(44) 


(45) 


For the baryon number density fluctuation Sn and the rapidity component of the fluid velocity SvP one can use the 
same expansion. For the fluid velocity components Su r and Su 4 we write instead 


Su r (r, r,(f>,p) = 2= [Su (t, r, </>, 77 ) + Su + (t, r, </>, 77 )] , 

fe 0 (r,r,^,77) = —= [Su~{T,r,<j>,p) - 5u + (t, r, 0, 77 )] 
r\J 2 


(46) 


with Su + * (r,r, <f>,p) = Su (r,r,<p,p). We expand Su [r,r,(j>,p) and 5u + {t, r, <j>, p) similar to Eq. (43) but replace 
Jm(kr) by J m _i(fcr) and J m +i(kr), respectively. The reality constraint becomes 


Su + *(T,k,m,q) = (—l) m+1 dM (■ T,k,—m,—q ). 


(47) 


In terms of the Bessel-Fourier transformed variables one can easily perform the spatial derivatives in Eqs. (38) - (|4 
To that end it is useful to use the relations 

Jm(kr) — — [</ m _i(fcr) T T m ^_i(/cr)], 
r 2 

d k 

— J m (fcr) = - [J m Jkr) - J m+1 {kr)}. 


(48) 
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FIG. 3. (Color online) Wavenumber k! l r "' > = z^ / R as a function of the discrete radial wavenumber l and for different values of 
the azimuthal wavenumber m — 1 (lowest curve) tom = 4 (uppermost curve). These values arise for the boundary condition 
Se = 0 at r = R and we choose R = 10 fm for definiteness. The plot shows that k increases with increasing values of both l and 
m, corresponding to finer spatial resolution. 


The Bessel expansion we use in Eqs. (43) contains an integral over all (positive) values of k. This expansion, also 


known as the Hankel transformation, is appropriate for functions on the open interval r € (0, oo). More realistically, 
the energy distribution in a heavy ion collision is non-zero only on a compact interval (0, R) with some radius R that 
depends on time during the expansion of the fireball and it is of the order of R ~ 10 fm. On such a compact interval 
the Bessel expansion becomes discrete, in the sense that the integral over k is replaced by a sum over a discrete subset. 
For example, the boundary condition Se = 0 at r = R leads to the values k\ rn> = Z( /R where the z < l r " > are the Z’th 
zero crossings of the Bessel function J m (z). To relatively good approximation zj m ^ is linear in m (for fixed l). In Fig. 
[ 3 ] we illustrate the resulting values for k\ m ^ as a function of the discrete radial wavenumber l and for different values 
of to. More generally, one might use an expansion based on J m ^( m )p(r)) where p(r) is a monotonous function into 
the interval (0,1) and a particularly useful choice for p(r) is discussed in Appendix A of Ref. [KTj . 

The evolution equation for the perturbation in energy density, Eq. ( |38[ ) becomes in Bessel-Fourier space (all per¬ 
turbation functions have now the argument (r, k , to, q) that we suppress for better readability) 


d T Se + 

+ 

+ 


- + -( a I 

T T \ 06 

1 (dp 
t \dn 


1 fdC 

de 


1 

t 2 V dn 


4 f dr/ 
3t 2 \ de 


4 fch 7 
3t 2 I dn 


Se 


Sn 


(49) 


2- 4 

e+p- - C+ z~y 

T 6T 


( Su + — Su ) + iqSvJ^J — —y iqSu 11 = 0. 


Similarly, the evolution equation for the perturbation in baryon number density becomes 


d T Sn -\—Sn + 

T 


nT 


e+p 


e + p 


nT n 2 


nT n2 


e + p 




AAin iws 


—j= (5u + — Su ) + iq Su v 


+ p2 1 6e 


(50) 




Let us now turn to the perturbations in the fluid velocity. Equations (40) and (41) lead to the following equations for 
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Su + and Su in Bessel-Fourier space 

2 


e+p- -C + —?7 

r St 


d T Su d 


=F 


dp 

Ik 

k k 2 


1 

T 

l^ 2 


dQ 

de 


2 

3r 


d T p-9 r C + 

T 


dr 1 

Ik 


1 - 2 4 

-^C + ^<9 r ?7+ 

r z Sr 


k , 

j 71 * t 


9p 

dn 


3t 2 

<K 

dn 


Su 1 


2 

3r 


Su^ — 

kfc 2 + k * 2 

Su T ± i 

\kq 

pkq 


.2 6 


_V2 



c?77 

9n 


<Su’ ! = 0, 


A*. 

J V2 


(51) 


and for the rapidity component we find from Eq. (421 


e+ p — —C — xk ) d T Su v + 

t St 


1 


1 — 


d T p+ -(e+p) d T (+ — C- 7T d rV- 


dp 

de 


-[C+i, 


1 

T 

ig/c 


t 2 \[2 


2 

3t 


977 

de 


*<7 


<5e + 


dp 

dn 


St 

dC 

dn 


3 r 2 




2 

St 


dr] 

dn 


%6n 

T 2 


(<5m + — 6u ) + 


C+ oV + ^ 


fe" = 0. 


(52) 


Note that eqs. (491 - (52) are now coupled ordinary differential equations. All spatial derivatives have become algebraic 
and one can directly integrate for the time dependent perturbations Se(T,k,m,q) etc. To construct such a solution 
one needs as an input the background or Bjorken solution for T(r) and p(r) as well as the relations that express 
all other thermodynamic densities (e, p, n etc.), transport coefficients (£, 77 , R) and derivatives (( dp/de ) n , ( dp/dn ) e , 
(i dC,/de) n etc.) in terms of the independent thermodynamic variables T and 77 . 

Let us first discuss some limiting cases of Eqs. (491 - (52) with extended symmetries. 


A. Statistical baryon number conjugation symmetry 


If the baryon number density vanishes in the background solution, i. e. n = /7 = 0, one has an extended symmetry 
namely baryon-anti-baryon or baryon number conjugation symmetry corresponding to , n —> —n. Odd derivatives 


such as ( dp/dn) t or ( dp/dn) t have to vanish and one finds that 5n decouples from the equations for Se in Eq. (49) 
and the perturbations of fluid velocity in Eqs. (511 and (52). However, this does not imply that Sn has to vanish as 


well. Locally and event-by-event one may have a non-zero baryon number density. The evolution equation for this 
perturbation is obtained from Eq. (50) as 


d T Sn -\—Sn + n 

T 


nT 


-P 




(53) 


The second term on the left hand side accounts simply for the dilution due to the longitudinal expansion while the 
third term is a diffusion term due to heat conductivity. Note that R is expected to be singular in the limit n —0 in 
such a way that the combination of terms that multiplies (k 2 + ^ 2 ) Sn remains finite [7S]- Therefore, the diffusion 
term indeed plays a role for the evolution of perturbations Sn. 

Equation (531 can be directly integrated and its solution reads as 


Sn(T,k,m,q) = (J^j exp [-fc 2 /i(r, r 0 ) - q 2 I 2 (T,T 0 )] Sn(T 0 ,k,m,q), 


(54) 


where the integrals 


h(r,T 0 ) 

h{r,T 0 ) 




d(p/T) 

dn 




( d(p/T) \ 

V dn ): 


(55) 


depend on the heat conductivity and thermodynamic quantities on the background Bjorken solution. While the 
integral I\ is typically dominated by late times r (for example for the ideal thermodynamic equation of state (14), 
heat conductivity of the form (22) and Bjorken expansion as in Eq. ©), the integral I 2 is dominated by early times 
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r ss To- Moreover, for fast thermalization —> 0 one has formally I 2 —> 00 such that in reality it might be rather 

large. Modes with q 7 0 are therefore strongly damped by dissipative effects of heat conductivity. 

The evolution equations for the perturbations in energy density Se and fluid velocity are independent of Sn. Their 
solution has already been discussed in a similar setup in Ref. @B]. 


B. Exact Bjorken boost symmetry 


The evolution equations for perturbations (491 - (521 simplify also in a situation where Bjorken boost invariance is 
realized as an exact symmetry instead of only on a statistical level. In that case one has Su v = 0 and the perturbations 
Se, Sn etc. vanish except for q = 0. Equation ((491) becomes 


d T Se - 


t t \oe 
1 (dp 
dn 


1 

T 2 

ac 

dn 


di] 

3 t 2 \Y 


Se 


drj 

3t 2 V dn 


Sn 


(56) 


2- 4 

e+p -C + 


t 3 r J y/2 


( 5u + — 5u ) = 0, 


and similarly Eq. (50) becomes 


d T Sn H —Sn + 

T 


nT n2 


*(£) 


(Su + — Sn ) 

y/2 V ' 


+ K 


nT n2 


m*- 


nT n2 


e + p 




(57) 


One observes that (56) and (57) depend on Su + and 8u only via the combination ( 5u + — 5u ) /y/2, for which one 
obtains from Eq. (51), 


(e+p- ^C+ (Su + -Su ) 


+ 


+ 


d T p - -d T C + + +Z d rV + 


3 r 


3r 2 ''J ^2 


-=(Su+-Su ) 


dp 

Ye 


1 (d( 

Ye 


drj 

3 t\Y 


kSe — 


dp 

dn 


1 (d C 

dn 


drj 
3 r V dn 


(58) 


k5n 


(k 2 + - rjk* 


y/2 


(■ Su + — Su ) = 0. 


Equations (561 - (58) together with the information about background quantities form a closed system that describes 
the analog of sound propagation and baryon number diffusion in the transverse plane of a longitudinally expanding 
fireball. The orthogonal combination of fluid velocity perturbations Su + + 5u~ is a shear mode with purely dissipative 
behavior (equation not shown). 

It is interesting to compare these equations to the ones that govern perturbations in a static medium. In that 
case all terms that involve explicit factors 1/r or derivatives of background quantities with respect to r vanish. For 
example, the analog of Eq. (56) is 


d T Se + (e + p) (<5m + — Su ) = 0, 
v 2 


(59) 


while the analog of Eq. (57) is 


d T Sn + n—j=(5u + — Su ) 


nT 




nT n2 


m/—■ 


(60) 
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/j/s = 1/(4 tt) nls = 10/(4tt) 





/J/'s - 1 /(] tt) /j/s = 10/(47t) 




FIG. 4. (Color online) Evolution of perturbations in energy density, baryon number density and fluid velocity with exact 
Bjorken boost symmetry (q = 0, Su v = 0) for different values of the azimuthal wavenumber m and radial wavenumber l. For 
R = 10 fm/c one has fcj 1 - 1 = 0.38fm -1 (black curves), k^ = 0.84fm -1 (red curves) and k^ = 1.30fm” 1 (blue curves). We 
compare two different values of the ratio of shear viscosity to entropy density r//s = 1/(4 tt) (left column) and (b) 77 /s = 10 /( 47 r) 
(right column). Heat conductivity is related to this by eq. ( |25| ). We use To = 0.5 GeV, /iq =0.05 GeV, to= 1 fm/c, r/=10 fm/c 
and for the initial values of the hydrodynamic fluctuations we choose Se(r 0 ) / 0, Sn(ro) = Su + (ro) = 8u~ (ro)=0. We denote 
A~ = u + — u~ (thus, Aq = Ug — Uq =0). See text for further details. 


and the analog of Eq. (58) is 


(e+p)d T ^= (6u+-6u k5e ~{j^) kSn 


- 4 ; 

C+3^7 


A: 2 — (<5m + — 5u ) = 0. 


( 61 ) 
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10 


/ 7 /s - 1/(4 tt) 



ql s= ! 0/(1 tt) 



FIG. 5. (Color online) Same as Fig. HI but for different initial values of the fluid perturbations: Sn(ro) ^ 0, <5e(ro) = <5m + (to) 
<5m - (to)=0 (thus, A,) - = Ug — Uq = 0). See text for further details. 


The set of equations (59) - to (61) describes sound propagation in the presence of dissipation due to shear viscosity, 
bulk viscosity and heat conductivity. We observe that at least some of the additional terms in Eq. (56) compared 
to (59) have the effect of an additional damping, in particular the square bracket in the first line of Eq. (56) is 
expected to be positive in the regime where fluid dynamics is applicable. Similarly, the leading additional term in 
Eq. (57) compared to Eq. (60) is the term ^Sn that has a damping effect, as well. The situation is less clear for the 
additional terms in Eq. (58) compared to Eq. (61), in particular the second line in Eq. (58) might actually conteract 
damping because d T p is negative. However, at least for larger values of the wavenumber k and non-zero viscosities 
the dissipative damping term in the last line of Eq. (58) is dominating. 

Equations (56) - ( |58| ) simplify further if the background is symmetric under baryon number conjugation as discussed 
in Sect. IV A In that case the perturbation in baryon number density 6n decouples from Eqs. (56) and (58) and is 
described by Eq. (53) (with q = 0). Nevertheless, the remaining equations for 6e and (Su + — Su ~) Jy/2 remain coupled 
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rjls = 1/(4tt) n/s= |o/( 47t ) 




ql s = 1/(4 tt) 




rjls - 10/(4tt) 



4 6 

r [fm/c] 


/ 7 /s = 10/(4 tt) 



FIG. 6. (Color online) Same as Fig. [4] but for different initial values of the fluid perturbations: <5u + (to)= 0.4, 8u (to)= 0.3 
(thus Ag = Uq — Ug = 0.1 and Aq" = Ug + Ug = 0.7), Se(ro) = 8n(ro) =0. See text for further details. 


and have to be integrated numerically for a given background solution and wavenumber k. This has already been 
discussed in Ref. [35] . 

In Figs. [4j [5] and [b] we show numerical solutions of the evolution equations (56 1 - (58) for the ideal EOS (14). For 
the background fields we employ the scaling solution ( |30| . We compare the numerical results for different initial 
conditions and two different values of the ratio of shear viscosity to entropy density and assume £ = 0 for simplicity. 
More precisely, the left columns of Figs. @1 and [ 6 ] correspond to 77 /s = 1/(47 t), the right columns to 77 /s = 10/(47 t). 
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/j/s = !/(4 tt) 5c () * 0 



///s = l/( 47 r) <5e„* 0 



qls = 1/(4tt) <5n 0 * 0 



/j/s = 1/(4 tt) <5A 0 * 0 



/j/s = 1/(4 tt) <5n 0 * 0 



/ 7 /s = l/(4rr) <5A(|* 0 



FIG. 7. (Color online) Amplitude of the perturbations at 77 = 10 fm/c in units of the initial weight at r = 1 fm/c as a 
function of the k- wave number. We choose r//s = l/(47r), To = 0.5 GeV and fig = 0.05 GeV. In each panel we choose a 
non-vanishing value initially for one of fluctuating fields while the remaining ones are set to zero.The top, middle and bottom 
panel corresponds to different initial conditions 5t 0 yf 0, Sno ^ 0 and A ( 7 = Uq — Uq yf 0, respectively. 


In all cases, the heat conductivity is taken to be related to the shear viscosity by (251. We also compare different 
values of the radial wavenumber k = k!' l ' n> = z\ m ^/R. We choose R = 10 fm/c which corresponds to k\ r> = 0.38 fm” 1 


(black curves), ki 2 ' 1 = 0.84fm 1 


v 2 — u.uimi (red curves) and k-f’ = 1.30fm” 1 (blue curves). In all cases, the modes with larger 
k are damped more quickly as expected. In order to simplify the notation we use the abbreviation = 5u + — Su~ 
in Figs. [4j [5] and [6] 

In Fig. [4] we have chosen initial conditions with non-vanishing perturbations in energy density <5e(r) = <5eo while 
the perturbations in baryon number density Sn and fluid velocity <5u + , <5w” vanish initially. The pressure gradients 
associated with <5e induce sound waves with the typical oscillating behavior between Se and Su + — Su ~, modified by 
the longitudinal expansion and viscous damping. As expected, the oscillation frequency is larger for larger radial 
wave numbers k. The perturbation in energy density 6e induces also a small perturbation in baryon number density 
Sn at times r > tq. This is due to the linear mixing between the different fluctuating fields (<5e, 8u + and 5u ~) for 
non-vanishing background baryon chemical potential (we choose /io = 0.05 GeV). For p> = n = 0, the evolution 
equation for Sn would decouple from the other fluctuating fields as we discussed in the previous section. Because we 
solve linearized equations for the perturbations, the solution scales linearly with the initial value Seg. 
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In Fig. [5]we initialize with non-vanishing perturbation in the baryon number density Sn(ro) = Sno but set de(ro) = 
5u + (tq) = Su~ (tq) = 0. The mode excited in this way has essentially diffusive behavior. This is most clearly seen 
in the intermediate panel which shows the temporal evolution of Sn/Sno- There are no oscillations seen but simply 
a decay in amplitude which is faster for large values of k. This decay is mainly a consequence of heat conductivity 
(or equivalently, baryon number diffusion). In addition to the baryon number density perturbation, also a (small) 
perturbation in Se and Su + — 8u~ is excited for r > tq. This is again a consequence of the non-vanishing baryon 
number density in the background. The behavior of these perturbations is oscillatory, i.e., of sound type. 

In Fig. 6 ]we choose initial conditions with Aq = 8uq — Suq 7 ^ 0 while the perturbations Se and Sn vanish initially. 
This results again in sound propagation of the typical oscillating type. In Fig. [b] we also show the behavior of 
perturbations in the orthogonal combination A + = Su + + Su~ which is a shear mode whose decay rate is determined 
by shear viscosity 77 . The shear viscosity dependence of the decay rate for this particular shear mode can be obtained 
directly from the corresponding evolution equation (51). 

In Fig. [d we show the final amplitude of the perturbations in energy density (left column) and particle density 
(right panel) at t/ = 10 frn/c as a function of the k- wave number in units of the initial weight at time tq=1 fm/c 
for 77 /s = 1/(47 t) and different initial conditions of the perturbations of the fluctuating fields. This plot shows that 
some modes of the initial perturbations characterized by the k- wave number indeed survive the entire evolution of the 
system and at the same time, it also indicates the distribution of of the surviving modes at the scales of time relevant 
for the freeze-out surface. 0 The uppermost panel corresponds to the non-zero value for the initial perturbation in 
energy density de(r 0 ) = <5eo while the remaining fluctuating fields, dn 0 , Su^ and 6uq , vanish exactly. The middle 
panel corresponds to the case where Sn(ro) = Sno and Se 0 = 6uq = Suq = 0. The bottom panel corresponds to 
Aq = Suq — Suq 7 ^ 0 and <5eo = Sno = 0. For the sound wave type initial conditions (deo 7 ^ 0 or A(j" 7 ^ 0) the size 
of the amplitudes at ry present a damped oscillatory behavior while for the case when Sno 7 ^ 0 the fluctuation of 
the Se and A - present an oscillatory behavior while Sn shows a exponential type decay which is typical to diffusive 
processes. In all the cases we observe that essentially none of the modes survive for values of k > 2 fin - > 

We conclude this subsection by emphasizing again the observation that perturbations in baryon number density 
have a diffusive time evolution with a dissipation rate determined by the heat conductivity. For typical values 
corresponding to strong coupling behavior, the damping is rather strong but at least the modes with the smallest radial 
and azimuthal wave-numbers (small values of m and l) are not dissipated completely and could have experimentally 
observables consequences. 


C. Exact transverse translation and rotation symmetry 


One can also consider a situation with exact symmetry under translations and rotations in the transverse plane. In 
that case only perturbations with k = 0 are possible and the fluid velocities in transverse directions have to vanish, 
Su + = Su~ = 0. Again Eqs. (49)-(52) simplify substantially, albeit not to a point where they can be integrated 
directly. Specifically, Eq. (49) becomes 
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and the evolution equation for the perturbation in baryon number density (501 becomes 
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(62) 


(63) 


7 This can be understood directly when taking the Fourier transform of the fluctuating fields at r = Tf, e.g., Eq. ( |43| ) The distribution of 
the fluctuating field as a function of its k- wave number determines the distribution of this field in coordinate space and determines the 
location of the maximum allowed correlation length in coordinate space. 
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FIG. 8. (Color online) Evolution of perturbations in energy density, baryon number density and fluid velocity with exact 
transverse translation and rotation symmetry (k = 0, Su + = Su~ = 0) for different values of the rapidity wavenumber q: q = 1 
(black line), q = 3 (red line) and q = 5 (blue line). We compare two different values of the ratio of shear viscosit y to entropy 
density 77 /s = l/(47r) (left column) and 77 /s = 10/(47 t) (right column). Heat conductivity is parametrized by Eq. ( |25| ) . We use 
To = 0.5 GeV, 710 =0.05 GeV, to= 1 fm/c, r/ = 10 fm/c and for the initial values of the hydrodynamic fluctuations we choose 
8t{r 0 ) / 0, Sn(ro) = 8u 71 (to) = 0. See text for further details. 


Finally, the evolution equation for the rapidity component of the fluid velocity (52) becomes 
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FIG. 9. (Color online) Same as Fig. [8] but for different initial values of the fluid perturbations: <5n(ro) ^ 0, <5e(ro) = 5 u v ( to ) = 0. 
See text for further details. 


These equations simplify further in a situation with vanishing baryon number density, the numerical solution for 
this situation has already been discussed in Ref. [36]. The solution of the fluctuating fields is in general complex but 
subject to the reality constraints 5e*(r,q) = Se(r, —q) and similar for 5n and Su 71 . 

In order to gain some qualitative insights let us consider a simple equation of state e = 3 p while setting n = jl = 0 
and neglecting the effects of viscosities where they are sub-leading compared to other background terms. One can 
then derive for the variable S = 5e/e the equation 


d 2 T S + 


5 

3 r 


c + 477/3 \ 

e+p ) 


drS+^5 = 0. 


(65) 


This equation describes sound propagation in the longitudinal direction on top of the expanding Bjorken background 
solution. Both the expansion and the viscosities have a damping effect as can be read of from the term ~ d T S. The 
last term in Eq. (65) is due to pressure gradients and the actual driving term of sound propagation. It is somewhat 
different than in other situations because of the time dependence ~ 1/t 2 . 
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/j/s=1/(4tt) /7/s = 10/(4 tt) 




/7/s = 1/(47/) /j/s = 10/(4jt) 



T [fm/c] t [fm/c] 


FIG. 10. (Color online) Same as Fig. [ 8 ]but for different initial values of the fluid perturbations: Su v (ro) 7 ^ 0 , <5e(ro) = 5ti(t 0 ) = 0 . 
See text for further details. 


More general, the set of equations (62), 


(631 and (64) describe also baryon number density waves and diffusion in 

111 and flO 


the longitudinal direction. We show numerical solutions to the evolution equations (621 - (64) in Figs. | 8 _ 
for different initial conditions. As we proceed in Sect. |IV B| we compare two different values of the ration of shear 
viscosity to entropy density 77 /s = 1 /( 47 t) (left panel) and 77 /s = 10 /( 47 t) (right panel). For the background fields we 
use again the scaling solution (30). 


For Fig. [ 8 ]we choose only Se to be non-zero initially. Compared with the behavior of the transverse sound waves 
or sound waves in a static medium discussed in the previous section, the evolution of the resulting longitudinal sound 
waves is completely different. In particular, no proper oscillations are visible during the entire temporal evolution. 
Rather, one observes a decay in amplitude, in particular at early times. This effect of the longitudinal expansion is 
particularly strong for large values of q. At later times the damping actually weakens to the extent that amplitudes 
remain non-zero at the final time. Interestingly, the influence of viscosity on the time evolution of longitudinal 
perturbations is relatively weak. Some quantitative differences are of course visible between the left and right panels 
of Fig. [ 8 ] but qualitatively, the evolution is surprisingly similar. 
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qls = 1/(47 t) 6e 0 * 0 



qls = 1/(4tt) <5e 0 * 0 



qls = 1/(4tt) <5n 0 * 0 



rjls = 1/(4tt) <5n 0 * 0 





q q 


FIG. 11. (Color Online) Amplitude of the perturbations at r/ = 10 fm/c in units of the initial weight at r = 1 fm/c as a 
function of the g-wave number. We choose q/s = l/(47r), To = 0.5 GeV and Ho = 0.05 GeV. In each panel we choose a 
non-vanishing value initially for one of the fluctuating fields while the remaining ones are set to zero. The top, middle and 
bottom panels correspond to different initial conditions 5e o ^ 0, 5n o ^ 0 and Mq ^ 0, respectively. 


Figure [9] was obtained by selecting only Sn ^ 0 at r 0 . Again, we do not observe any proper oscillations of the 
fluctuating fields along the longitudinal direction. In this case the viscosity and heat conductivity have a somewhat 
larger effect. The amplitude of the fluctuating fields gets damped as one increases the value for the shear viscosity 
and heat conductivity (according to Eq. (25)). 

In Figure [To] we choose Suq ^ 0. As in the previous two situations, there is no proper oscillation visible for the 
time interval shown. 

Finally, Fig. [IT] shows the final amplitude of perturbations at r/ = 10 fm/c in units of the initial weight at time 
r 0 as a function of the longitudinal q- wavenumber. Figure [TT] is obtained by choosing a non-vanishing value initially 
for one particular fluc tuat ing field while the remaining fluctuating fields are initially set to zero. The top, middle and 


bottom panels of Fig. 11 corresponds to deg ^ 0, Sng ^ 0 and Su 0 ^ 0, respectively. We observe that the amplitude 


of the fluctuating modes goes asymptotically to zero for q > 25 which corresponds to a small window in the rapidity 
variable (i.e., Ar) ~ (Aq)” 1 ). We expect that modes with intermediate and large q would be damped stronger for 
earlier initialization time To- 

Finally, in a situation where Bjorken boost invariance as well as translations and rotations in the transverse plane 
are realized exactly, i. e., 5u v = 6u + = Su~ = k = q = 0, Eqs. (49) and (50) reduce simply to a linearized version of 
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the Bjorken expansion in Eq. (26) as it has to be. 


V. THE TWO POINT CORRELATION FUNCTION OF BARYONIC PARTICLES 


In this section we discuss the possibility to access the information about perturbations in the baryon number density 
experimentally by measuring a correlation function of the net number of baryons (baryons minus anti-baryons) as 
a function of the rapidity and azimuthal angle. We concentrate for simplicity on the case of vanishing background 
baryon number density as discussed in Sect. |IV A| 

Perturbations in baryon number density in position space as described by Eq. ( [54] ) are not directly accessible to 
experiments. However, a fluctuating baryon number density and chemical potential on the kinetic freeze-out surface 
has an influence on the distribution of particles with non-zero baryon number in momentum space. This concerns 
in particular protons but also resonances with non-vanishing baryon number. Similar as for flow observables, there 
is a direct link between different harmonics in azimuthal angle and rapidity in the fluid dynamic description and 
the corresponding harmonics in the momentum space particle distribution. Thus, we can partly access the physical 
information contained in Eq. ( |54[ ). As an example, we consider a connected two-point correlation function of the typ^] 

^Baryon(^l 4 > 2:Vl V 2 ) = (^Baryons(01) Vl )^Baryons (02 1 V2 ))c> (66) 

which measures correlations of baryonic particles (i.e. the number of baryons minus anti-baryons) as a function of 
the difference between (particle momentum) azimuthal angles 0i — 02 and (particle momentum) rapidities rj\ — 772 - 
In Eq. (66), nBaryons(0, v) is the number of baryons minus anti-baryons as found in the detector in a particular bin in 
azimuthal angle 0 and rapidity 77^] We also introduce the Fourier representation 


CBaryon(01 ^2? ^7l ^2) — ^ ^Baryon Q) 1 


,im((f > 1 — tf>2)+iq(ri 1 -ri2) 


(67) 


The correlation function in Eqs. (66) and (|67|) is determined by a combination of initial conditions (set at the time 


where a fluid dynamic description becomes valid) and response functions that describe how baryon number density 
perturbations propagate in the fluid dynamic regime and how they influence the particle distributions at freeze-out. 

In the following we discuss both parts in a bit more detail. First, the initial state after a heavy ion collision (and after 
the early non-equilibrium dynamics) at the time To when a fluid dynamic description becomes valid, is characterized 
by a fluctuating baryon number density (fn(T 0 , r, 0,77) around some average or expectation value n(To,r). (The latter 
might be rather small at LHC and upper RHIC energies and we neglect it in the following.) For the fluctuating part 
we use a Bessel-Fourier decomposition 


Sn(T 0 ,r,(j),ri) = 


E E 

m =—oo l—l ' 


dq 


(m) 


(q) e i ^+ i ”>4(z i W p(r)) 


( 68 ) 


An event-by-event ensemble of initial conditions conditions for the baryon number density can be characterized in 
terms of the weights Sn^ m \q). For example, the two-mode correlation function is 


K mi V) H ( 2 m2) (^)) = 2 tt%i + q 2 )S, 


mi+m2,0 C'SnSndi.lo (9)- 


(69) 

We have assumed here that the ensemble of initial conditions is statistically symmetric under azimuthal rotations and 
longitudinal boosts leading to the factors S mi+m2 ,o and 2ir5(qi + q 2 ) on the right hand side of Eq. (69). 

For a single event with baryon number perturbation as in ( |68[ ) , the baryon number distribution in momentum space 
after kinetic freeze-out will be proportional to the weights 5rij" l> (q) within the linear response approximation. More 
specific, the Bjorken-boost and azimuthal rotation symmetries imply that one can write 


Am) 
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(?) = E Aryans-,{m)l(q)Sn[ m) (q), 


(70) 


The brackets (• • •) in Eq. < |66| ) denote an event average 
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9 There is a complication due to the fact that neutrons cannot be measured experimentally. Further studies are needed in order to quantify 
whether this presents a problem for observables as in ( |66[ ) and if so, how these can be overcome. Also, one should estimate possible 
contributions to Eq. ( |66| ) from sources other than fluid dynamics, such as resonance decays. 
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with linear baryon number response function S , B ar yons;(m)i(<z)- The object on the left hand side of Eq. (70) is the 
Bessel-Fourier weight of the (momentum space) distribution of the number of baryons minus anti-baryons. The 
correlation function on the right hand side of (|67|) can be written as 


CBaryon(w, q) ^ ' ‘S'Baryon;(m)Zi (?) ^Baryon^— m)l 2 ( l) ^SnSn',l\ ,l 2 (®) 
=1 


(m) 


(71) 


h,i 


For a more detailed discussion of the response function formalism briefly introduced here we refer to Ref. [FT . 

The linear response functions SBaryon ; (m)z(?) are i n particular also affected by heat conductivity. More specific, the 
analog of the factor exp(— k 2 Ii — q 2 !^) in a situation with realistic transverse dependence and radial flow leads to a 
suppression of modes with q 2 > 0 and large values of m and/or the radial wave number l. Qualitatively, one expects 
that the scale for the suppression in the transverse direction is set by the (time dependent) radius R of the fireball. 
Moreover, the Z’th zero crossings z^'"' 1 of the Bessel-functions J m {z) are for fixed l approximately linear in m (for the 
relevant values of to and l , with prefactor of order unity) so that one expects qualitatively 

CBaryon(m, q) ~ exp(-2 m 2 1[ - 2 q 2 I' 2 )C£=° yon (m, q), (72) 

where on the right hand side Cg=° yon (m, q) would be the corresponding correlation function in the (somewhat hy¬ 
pothetical) situation of vanishing heat conductivity and the dissipative attenuation terms can be roughly estimated 
as 
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(73) 


Going now back to the two-particle correlation function (66), the exponential suppression factor in Eq. (72) implies 


for large I 2 long range correlations with respect to the rapidity difference 771 — 772 , with a decay that is determined 
by the value of I' 2 (except if C , Baryon( m > q) has a very strong decay with q already) and a similar, although weaker, 
effect with respect to the azimuthal wavenumber to. In order to make our qualitative statements more precise, it is 
necessary to generalize the calculations described here to a more realistic background. A more realistic background 
would have a realistic transverse profile and expansion in addition to the longitudinal (boost-invariant) expansion. 
Moreover, one also has to perform more detailed studies of the initial conditions and kinetic freeze-out, that both 
affect C|=° yon (TO, 9 ). 


VI. CONCLUSIONS 

We have studied solutions of the fluid equations describing relativistic heavy ion collisions in the presence of a 
globally conserved quantum number (baryon number) using a background-fluctuating splitting. For the background 
we have assumed Bjorken boost and transverse translation and rotation invariance. This generalizes Bjorken’s original 
solution to non-vanishing baryon number density as well as shear and bulk viscosities. Heat conductivity does not 
play a role on the background equations since the diffusion current vanishes exactly due to the symmetries of the 
Bjorken flow. 

We derived evolution equations for the perturbations around this background solution. While the amplitude of 
these perturbations was assumed to be small, such that linearized equations could be used, the formalism allows us to 
treat perturbations with arbitrary dependence on the transverse coordinates and rapidity. Technically, this is done by 
employing a Bessel-Fourier expansion. The partial differential equations of relativistic fluid dynamics become ordinary 
differential equations for the different modes that are characterized by radial, azimuthal and rapidity wave numbers. 
The evolution of these perturbations is governed by the thermodynamic properties encoded in the equation of state 
p(T, n) as well as the transport properties (i.e., shear viscosity bulk viscosity £(T, /i) and heat conductivity 

tc(T,/ r) in the first order formalism we use). 

Generically, one finds that perturbations with large wave numbers are damped more quickly by the dissipative 
processes, as expected. The dissipation of different modes depends on time in a different way and, in particular, 
deviations from Bjorken boost symmetry show a fast damping at early times. In principle, it might be possible to use 
these dependencies to probe transport and thermodynamic properties at different times in the evolution history and 
therefore for different temperatures of the quark-gluon plasma produced in a heavy ion collision. 
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In order to make more quantitative statements, one must take a realistic transverse density profile and expansion 
into account, of course. This has been done for perturbations with exact Bjorken boost symmetry and vanishing 
baryon number in Refs. 03 M E33 HU- In the present paper we have concentrated mainly on the evolution of 
perturbations in baryon number density. They have diffusion-type evolution governed by the longitudinal expansion 
and heat conductivity. (In the Landau frame, heat conductivity can in fact be understood as baryon number diffusion.) 
There are characteristic differences in the dependencies on longitudinal and transverse wave numbers. More specific, 
baryon number perturbations are quickly “flattened out” in the longitudinal direction at early times. 

In principle, the information on baryon number perturbations is accessible experimentally via two-point (and higher 
order) correlation functions of particles with non-zero baryon number, as a function of the difference in azimuthal 
angles and rapidities. Based on the evolution equations for perturbations, we expect long-range correlations in rapidity 
(a “baryon number ridge”). For a more detailed theoretical picture one needs a better description of the local event-by¬ 
event fluctuations in baryon number density at the initial time when fluid dynamics becomes valid. Also, one should 
take a realistic transverse expansion into account and study the implications of baryon number perturbations at the 
kinetic freeze-out. (Formulas needed for this have already been derived in Ref. [50].) It would be very interesting to 
study net-baryon number correlations experimentally, as well as theoretically in more detail, and thereby constrain 
heat conductivity as another property of the quark-gluon plasma. 
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Appendix A: Thermodynamic relations in the grand canonical ensemble 


In this appendix we compile some thermodynamic relations in the grand canonical ensemble that we found useful 
in the context of relativistic fluid dynamics with a conserved charge. We start from the pressure p(T,/it), which is 
related to the thermodynamic potential of the grand canonical ensemble (the Landau potential) by p = —f 1/V. The 
differential of pressure is 


dp = sdT + ndp,. 


(Al) 


All thermodynamic quantities can be obtained from this and the Gibbs-Duhem relation e +p = Ts + /m, for example 


s = 




(A2) 


In the following we will sometimes drop the subscripts with the convention that pressure is evaluated as a function 
of T and p unless indicated otherwise. Also we find it useful to express all susceptibilities in terms of the pressure 
and its derivatives. This avoids ambiguities and realizes Maxwell’s relations automatically. For example, the energy 
density is obtained then as 


, rp d p dp 
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Its differential, as well as the one for density, are 
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These linear relations can be inverted to yield dT and dp in terms of de and dn, 
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Other useful quantities are the heat capacity densities 
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the isothermal and adiabatic compressibilities 
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the thermal expansion coefficient 
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and a modified sound velocity at fixed particle density 
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Both sound velocities agree for vanishing baryon number density, n = 0. Note that the usual relations 
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are fulfilled. Moreover, one has 
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For the evolution equations of linear perturbations as discussed in Sec. IV we need also 
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(A13) 


where, similarly to pressure p(T, //), the bulk viscosity £(T, p) and shear viscosity rj(T, p) are functions of T and p on 
the right hand side. 


Appendix B: Linearized relativistic fluid dynamics 


In this appendix we discuss a background-fluctuation splitting for the fluid dynamic equations 
Sec. |IV| We split the fluid dynamic fields into a background part and a perturbation according to 

u 11 = + Su M , 

e = e + Se, 
n = n + Sn, 

VCbl 1 ]k — Tibulk “t" d7Trj111k. 


as it is used in 


(Bl) 


and so on. The projector orthogonal to the fluid velocity is given by 

AP V = Kp v + SA^, (B2) 


with 6A >J ’ U = u^Su 12 + Su^u 12 . If one restricts to a linear treatment of perturbations, the equations of motion for the 
background are simply the full equations of motion ([5]). For the perturbations in energy and particle number density 
one obtains from Eq. ([5]) 
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In these equations one can see 5e and Sn as independent variables, to which other thermodynamic variables and the 
transport coefficients are related in the standard way, e.g., 


Sp 

d^Sp 



a " ie+ GI) a r 6 n + a. 



Se + d^ 
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From the constitutive relation of first order fluid dynamics in Eq. (p| one finds 
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Similarly, for the bulk viscous pressure in Eq. 0 one finds 
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Finally, the perturbation of the diffusion current is obtained from eq. (|8]) as 
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Equations (B3) and (B4) also involve the following divergence of the shear stress perturbation 
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the derivative of the bulk viscous pressure perturbation 
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dpSO = V M V a <5w“, 

and finally the divergence of the perturbation in the diffusion current 
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+ V, 




Note that the expression in Eq. (B19) is contracted in Eq. ( |B18 ) wi th th e projector A“^. In many circumstances the 
background field changes only in the direction of such that Eq. ( |B19 ) simplifies substantially. 

Note that the formulas compiled in this appendix allow us to obtain for a given background solution and ther¬ 
modynamic equation of state linear evolution equations for the perturbations around this background solution. The 
independent variables of these linearized equations are the three independent components of Su ^ (one constraint is 
given by the condition il^Su 11 = 0) as well as 6e and Sn. In the first order formalism of relativistic fluid dynamics, 
the equations for the perturbations are of parabolic type while they are expected to become of elliptic type when 
relaxation time terms are kept. 
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